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Abstract 

In this paper we present a Path Integral Monte Carlo (PIMC) simulation 
of the orthorhombic phase of crystalline polyethylene, using an explicit atom 
force field with unconstrained bond lengths and angles. This work represents 
a quantum extension of our recent classical simulation [|^. It is aimed both 
at exploring the applicability of the PIMC method on such polymer crystal 
systems, as well as on a detailed assessment of the importance of quantum 
effects on different quantities. We used the NpT ensemble and simulated the 
system at zero pressure in the temperature range 25 - 300 K, using Trotter 
numbers between 12 and 144. In order to investigate finite-size effects, we used 
chains of two different lengths, C12 and C24, corresponding to the total number 
of atoms in the super-cell being 432 and 864, respectively. We show here 
the results for structural parameters, like the orthorhombic lattice constants 
a,b,c, and also fluctuations of internal parameters of the chains, such as 
bond lengths and bond and torsional angles. We have also determined the 
internal energy and diagonal elastic constants cii,C22 and C33. We discuss 
the temperature dependence of the measured quantities and compare to that 
obtained from the classical simulation. For some quantities, we discuss the 
way they are related to the torsional angle fluctuation. In case of the lattice 
parameters we compare our results to those obtained from other theoretical 
approaches as well as to some available experimental data. In order to study 
isotope effects, we simulated also a deuterated polyethylene crystal at low 
temperature. We also suggest possible ways of extending this study and 
present some general considerations concerning modeling of polymer crystals. 
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I. INTRODUCTION 



Polymer crystals are well known to be intrinsically difficult to prepare in a highly crys- 
talline state, which in turn hinders the possibilities of their experimental characterization. 
As a consequence, computer simulation appears to be a convenient tool to study their prop- 
erties. For crystalline PE, which represents the simplest and thus paradigmatic case, it has 
been recently shown |l|] that classical constant pressure MC is a well applicable simulation 
method, provided a good quality force field is available. It allows to calculate the whole 
variety of static local and collective quantities, including properties of major practical and 
technological importance, like thermal expansion and elastic constants. On the other hand, 
recent work on the same system, using quasi-harmonic or self-consistent quasi-harmonic ap- 
proximation 0-0, has clearly pointed to a quantitative as well as a qualitative inadequacy 
of the classical treatment at low temperatures, where quantum effects cannot be neglected 
anymore. 

Generally, quantum effects are known to be important in lattice dynamics of solids in the 
low-temperature region, when classical thermal fluctuations become comparable or smaller 
than the amplitude of the quantum zero-point motion. Under particular circumstances, 
when two different crystal structures have classically very close energies, and the system is 
close to a structural phase transition, quantum fluctuations can play a decisive role. In case 
of solid nitrogen 0, they are responsible for a strong isotope effect on the low-temperature 
a — 7 structural phase transition as a function of pressure. If the classical energy difference 
is small enough, quantum effects can even suppress the transition altogether, and stabilize 
the disordered phase, as it happens in the case of quantum paraelectrics SrTiOs or KTaOa, 
where the ferroelectric long-range order is only incipient down to T = 0. Even though a 
PE crystal does not represent such a dramatic case, at low temperatures it is not possible 
to account even qualitatively for the temperature dependence of quantities like thermal 
expansion coefficients and elastic constants without quantum effects being duly taken into 
account. 

A natural extension of the classical MC method in order to include quantum effects at 
finite temperature is the Path Integral Monte Carlo scheme [|^. Recently, the NpT version 
of this method has been applied to study quantum effects in crystals at low temperatures, 
in particular in solid rare gas systems and silicon crystal In case of polymer crystals, 
the distinguishing feature is an extreme anisotropy, closely related to the existence of many 
energy scales, ranging from soft intermolecular (non-bonded) interactions to stiff intramolec- 
ular (bonded) interactions. This feature presents a problem already at the classical level, 
as discussed in |1T0| , pTl ,p]| , and requires an introduction of special global moves in the sam- 
pling algorithm. In the quantum case, we should moreover expect very different convergence 
properties of different physical quantities as a function of Trotter number, depending on the 
typical energy scale with which a given quantity is associated. 

The aim of the paper is basically twofold. On the one hand, we wanted to explore the 
applicability of constant pressure PIMC method to a PE crystal, and determine the region 
of temperatures where the use of the method is practical. On the other hand, since the 
PIMC scheme is capable of providing essentially exact results, it can also be used to assess 
the range of validity of approximate analytical methods, like, e.g., quantum quasi-harmonic 
approximation . This is particularly important for the study of intrinsically anharmonic 
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phenomena, like, e.g., lattice thermal expansion. 

The paper is organized as follows. In section ||, we briefly describe the PIMC simulation 
method used, without addressing the force fleld and its implementation, since these issues 
have already been discussed in detail in |l|]. In section |T|, we present and discuss the results, 
paying particular attention to comparison of the quantum results to the classical ones. In 
the final section we draw some conclusions and suggest some possible ways of extending 
this study. We also make several remarks concerning general issues related to modeling of 
polymer crystals. For completeness, in Appendix A we present the full form of the force 
field we have used together with numerical values of the parameters. In Appendix B, we 
present some considerations on the relation of correlation functions of torsional fiuctuations 
to the contraction of the crystal along the chain axis. 



II. PIMC SIMULATION METHOD 

In this section, we will describe only those features of the simulation method which are 
specific for the constant pressure PIMC scheme. The implementation of the SLKB force field 
we used as well as many other features of the present algorithm are exactly identical 



to those of our classical simulation, which has been described in detail in Ref. [I]]. For 
convenience, however, in Appendix A we summarize the form and parameters of the force 
field. 

We have implemented the constant pressure PIMC scheme basically along the same lines 
as it was done for a cubic system in Ref. |^, the only difference being that in our case 
we had to use an anisotropic version of the NpT ensemble. We have used the primitive 
decomposition of the hamiltonian, resulting in the effective hamiltonian 

Hp{m = E (e (-t - -rf + ^^({-t})) , (1) 



fc=i 



where is the number of particles in the quantum system, rrii are their masses, P is the 
Trotter number, j3 = l/ksT is the inverse temperature and is the potential energy 

of the system. Such an effective hamiltonian represents a pseudo-classical system consisting 
of P copies (Trotter slices) of the original system, individual particles in neighboring Trot- 
ter slices being connected via harmonic "springs", and periodic boundary conditions being 
applied along the Trotter direction. This pseudo-classical system has now NP particles 
and can be simulated using the same constant pressure MC algorithm as in the classical 
case. The acceptance criterion for the volume moves was based on the Boltzmann factor 
(■5iS2S3)^^e~^^^, where Ep = Hp + PV0S1S2S3, p is the external pressure, Vq is the volume 
of the reference super-cell and si, S2, S3 are three independent scaling factors along the coor- 
dinate axes. In all simulations described in this paper, the external pressure was set to zero. 
The estimators for all configurational properties, diagonal in the coordinate representation, 
are straightforward analogues of their classical counterparts, while for the kinetic energy we 



used the virial estimator [T^. 



To sample the system, we have used three kinds of moves: classical moves, quantum 
moves and volume moves. Classical moves of two types, local moves of atoms or global 
moves of whole polymer chains, as described in Ref. [|lO|,|Tl],|l| , have now always been applied 
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to all particles or chains with a given number in all Trotter slices simultaneously. These 
moves sample the classical configurational phase space of the system. We note here that 
when performing a rotation of a given chain in all Trotter slices, the energies of the "harmonic 
springs" between the corresponding individual particles have to be recalculated explicitly, 
in contrast to pure translational moves of a particle in all Trotter slices, which preserve the 
energy of the "springs" . Quantum moves consisted of local translational moves of individual 
particles of the pseudo-classical system, which sample the quantum fluctuations around the 
classical paths. In the quantum moves, different maximum displacements have been used for 
C and H atoms, not only because of the different number of bonds but also because of the 
different mass of the atoms and resulting different stiffness of the "springs" . In volume moves, 
we performed a simultaneous rescaling of coordinates of all particles in all Trotter slices by 
the three scaling factors si, 52,53. One Monte Carlo step per site (MCS) thus consisted of 
an attempted quantum move on each particle of the pseudo-classical system, followed by 
a classical move attempted successively on all atoms or all chains (always simultaneously 
in all Trotter slices, as described above) and a volume move. Among the classical moves, 
30 % of global moves were used, like in the classical study For all kinds of moves, the 
displacements were chosen to yield an acceptance ratio of 20 - 30 %. 

We have simulated systems with C12 and C24 chains, consisting of 2 x 3 x 6 and 2 x 3 x 12 
unit cells, respectively (432 and 864 atoms), at temperatures 25, 50, 100, 150, 200 and 
300 K. At different temperatures, different numbers of values of the Trotter number were 
used. For the smaller system, we used at 25 K P = 144, at 50 K P = 54, 72, 144, at 
100 K P = 36, 54, 72, at 150 K P = 48, at 200 K P = 16, 24, 32 and finally at 300 K 
P = 12, 16, 24. The larger system was simulated only at 100 K with P = 72 and at 300 K 
with P = 24. We note here that the largest pseudo-classical systems simulated consisted 
of 432 X 144 = 864 x 72 = 62208 particles, which is quite a large number. As an initial 
configuration for a given temperature, we always used a pseudo-classical system consisting 
of P identical copies of an equilibrated classical configuration at the same temperature. 
This configuration was then equilibrated for several thousands MCS with the full PIMC 
algorithm, which corresponded to switching on the quantum fluctuations and allowing the 
system to find a new equilibrium. For illustration of the run length used for measurement, 
for the smaller system at 300 K and P = 24 we used about 280000 MCS. Roughly the 
same amount of CPU time was used for all data points, the number of Monte Carlo steps 
per site thus scaling inversely with number of particles NP of the pseudo-classical system. 
The whole run consisted of numbers of subbatches ranging from 5 to 57 and the batch 
subaverages were used to estimate the approximate error bars of the total averages. 

III. RESULTS AND DISCUSSION 

Before discussing in detail the results for various quantities, we comment briefly on 
convergence of the PIMC scheme as a function of the Trotter number P. For different 
quantities, we have found considerably different convergence, the best case being that of 
quantities like, e.g., lattice constants, which are mainly related to softer interactions (non- 
bonded interactions and torsional terms). In such cases, where the results obtained with 
different values of Trotter number P were identical within the statistical error, the quantum 
limit was practically reached and no extrapolation was necessary. A considerably slower 
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convergence is found for quantities like the energy, which depend crucially on fluctuations 
of degrees of freedom related to strong (bonded) interaction potentials. In these cases, an 
extrapolation to P — oo was performed in order to recover the true quantum values. We 
have used the standard formula [1141 



^P = ^oo + |^ + ^ + 0(-l), (2) 

which requires data for three different values of Trotter number P in order to find the 
extrapolated value A^o. 

In the discussion, we concentrate mainly on the comparison of quantities obtained from 
the quantum simulation to their classical counterparts, presented in We start with the 
local quantities, in particular fluctuations of the internal coordinates, for which the quantum 
effects are found to be most pronounced. 

In Figs. and ^, we show the temperature dependence of the average bond length 
fluctuation for C-C and C-H bonds, respectively. The distinguishing feature of quantum 
results for such fluctuations is their saturation at rather large values at low temperatures. 



instead of the classical vanishing {\J {{5rY) oc -\/T as T — > in the classical case). In 
particular in case of the C-H bond, there is a marked Trotter dependence of the results. By 
means of the above described extrapolation, we found at T = 50 K the values of 0.05 A for 



{{5rccY) and 0.079 A for y {{SrcuY), which represent about 3 % and 7 % of the respective 
equilibrium bond length. These values are representative of the ground state values, since 
both curves appear to be entirely flat up to room temperature, reflecting high frequencies of 
corresponding bond stretching phonon modes. In Figs. ^ and |[ average fluctuations of the 
bond angles Occc and 9hch are shown as a function of temperature. Trotter extrapolation 



at T = 50 K yields here the values of 3.46° and 8.44° for ^{{decccY) and ^{{deHcnY), 
respectively. While the former curve increases at room temperature by about 0.5° with 
respect to the ground state value, the latter one corresponding to the bond angle involving 
two hydrogen atoms is completely flat. 

In Fig. ^, we show the temperature dependence of the average torsional angle fluctuation 
i^cccc) froKi the trans minimum, which already in the classical case has been shown 
to play an important role in the physics of the system. The values at T = 50 K and T = 25 
K are already very close to each other, and therefore the T = 25 K value of 5.59° can be 
considered to be representative of the ground state. The characteristic temperature, below 
which the difference between the classical and quantum values starts to increase rapidly, can 
be estimated to be about 150 K. At room temperature, the difference is still about 0.5°. 

The internal energy per unit cell of the system is shown in Fig. The dependence on P 
is in this case particularly pronounced, since the dominant contribution at low temperatures 
comes from the hard degrees of freedom which require larger values of the Trotter number 
in order to converge. The extrapolation to P ^ oo is thus absolutely necessary here, and 
has been performed for all temperatures where data for three values of P are available, 
i.e. 50, 100, 200 and 300 K. The extrapolated values at 50 K and 100 K suggest that the 
energy in this region would still somewhat decrease by approaching zero temperature and 
the extrapolated value of 62.125 kcal/mol at 50 K can be considered as an upper estimate 
of the ground state total energy. By subtracting the classical ground state energy which is 
found to be equal to -7.39 kcal/mol per unit cell, we flnd a value of 17.38 kcal/mol per CH2 
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group, which agrees well with the zero-point energy of 17.598 kcal/mol, obtained in Ref. |T5 
within quasi-harmonic approximation for a different force field. 

Concerning the structural stability, for the smaller system we observed at T = 300 K an 
occasional rotation of a whole chain from the "herringbone" structure, just as in the classical 
case [|I|. No such rotation was observed for the larger system, however. 

Before passing to the discussion of the temperature dependence of the lattice constants a, 
b and c, we would like to make a general comment on the statistical error of lattice constants 
of crystals evaluated within constant pressure PIMC scheme. The statistical error a{{a)run) 
of the lattice constant {a)run averaged over a run consisting of N configurations is given by 

,, , , cr(fl) /„s 
(J{{a}run) = —r^ , (3) 



where a{a) is the intrinsic fluctuation of the quantity a, and s is the corresponding statistical 
inefficiency, expressing the effect of correlations between subsequent configurations of the 
Monte Carlo run |16|. In order to keep the Trotter error roughly constant at different 



temperatures, one usually keeps the product PT constant. For given system size and amount 
of CPU time, the number of configurations N is inversely proportional to the Trotter number 
P and thus directly proportional to temperature T. On the other hand, the fiuctuation 
a (a) of the lattice constant a is essentially a fiuctuation of the linear size of the system 
which is a purely classical fiuctuation, and obeys the equipartition theorem. Provided the 
elastic constants of the system do not vary too much, which should be well satisfied at low 
temperatures, the relation cT^(a) ~ T should hold. Combining the expressions together, we 
find 



(^{{a)run) (4) 



where the explicit dependence on T in numerator and denominator has cancelled. Of course, 
there is still an implicit dependence hidden in the factor s, which increases with decreasing 
temperature because of growing system size in the Trotter direction. Nevertheless, in case 
of such purely classical fiuctuation which vanishes as T — > 0, the situation is more favorable 
with respect to the case of a general fiuctuation which would instead tend to a finite zero- 
point quantum value. In our results for the lattice constants, this is illustrated by the fact 
that the error bars of points at lowest temperatures are not larger than those of points at 
higher temperatures. 

The temperature dependence of the lattice constant c is shown in Fig. |^. In both, 
classical and quantum lattice contraction with increasing temperature is observed. 

An interesting feature here is that at T ~ 50 K, the classical and quantum curves cross. At 
higher temperatures, the quantum result stays above the classical one while falling slightly 
below that for T < 50 K, where the quantum fiattening appears. This behavior suggests 
the presence of at least two distinct quantum effects. Since in the classical case |^ it turned 
out that there is a linear dependence between c and {(fcccc)^ have tried the same plot 
for our quantum data. Fig. ||, in order to separate different contributions. In the latter plot, 
the classical and quantum curves are roughly parallel to each other, which shows that the 
low-temperature flattening of c is a consequence of freezing of the thermal contribution to 
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the torsional fluctuations. With the exception of the lowest temperatures, the shift of the 
quantum curve with respect to the classical one appears to be temperature independent, 
being roughly equal to O.OOISA. This represents a zero-point expansion of the lattice along 
the c-direction arising from hard modes which are not significantly excited even at room 
temperature. While in the classical case all points are found to fall very well on a straight 
line in the whole temperature range in the quantum case a distinct upward bending of 
the curve is apparent at low temperatures. To understand its origin, one might try to use 
the full formula 
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derived in Ref. (see also Appendix B), which relates the contraction of the lattice con- 
stant c to correlation functions of two neighboring torsional angle fiuctuations 0o, 0i. The 



temperature dependence of the normalized correlation function 



(</'O0l) 

1W 



is shown in Fig. 0. In 



contrast to the classical curve, which is fiat in the whole temperature region, the quantum 
one is seen to increase in value at low temperatures where it becomes almost a factor of 
two larger than the classical one. In Fig. |10| we plot c vs. |((0o — In this plot, the 

low temperature upward bending from Fig. ^ has become less pronounced and the points 
now exhibit a clear linear dependence, which suggests that the bending has its origin in the 
quantum reinforcement of the correlation function -^^^Ij^ of neighboring torsions at low tem- 
peratures. The slope of 1.275 x 10~^A deg~^, however, turns out to be larger in magnitude 
by a factor of two with respect to the value of |cosin^ f (ifo)^ ~ ^-^^ ^ 10~^A deg~^ result- 
ing from Eq. (we took the values c = 2.53A and a = 180° — 110.75°). About the same 
discrepancy in the slope is found also for the classical data, and this has been pointed out 
already in Ref. [|T| (where, however, the correlation function between neighboring torsions 
was neglected and the discrepancy thus turned out to be just about a factor of 1.5). While 
in the classical case the thermal contraction could be modified due to contribution of modes 
other than torsions, in the quantum case such harder modes are mostly frozen even at room 
temperature, and do not contribute considerably. Thus in the quantum case the formula (^ 
should yield a better agreement with the simulation. As it turns out, however, it captures 
qualitatively correctly the basic role of torsional correlation functions in the thermal con- 
traction, but fails in the quantitative aspect. We discuss the origin of this discrepancy in 
detail in the Appendix B. Analogously to the classical case, no significant finite-size effects 



are seen on the lattice constant c. Comparing the quantum result to experimental data |T7 
in Fig. 1^, we see that apart from the constant offset, the agreement has improved at low 
temperatures, due to the quantum flattening. 



In Fig. 11, we show the temperature dependence of the lattice constant a. At tempera- 
tures below 50 K, the quantum curve appears to be entirely flat, and the difference between 
the quantum and the classical result is in this region as large as O.ISA, which represents a 2 
% effect. At all temperatures, the quantum curve lies above the classical one. It is interesting 
to note that the difference between the two curves persists up to room temperature, being at 
T = 300 K equal to 0.06A, which is still about a half of the zero-temperature value. For this 
lattice constant, a very good agreement between the simulation and experimental results 



1% is found, which proves that almost the whole low-temperature discrepancy between the 



classical simulation results and experiment is purely due to quantum effects. Similarly to 
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the classical case, no finite-size effect is seen on the quantum curve at T = 100 K, while a 
small one is seen at T = 300 K. 



Analogously to the case of the lattice constant c, in Fig. |12] we plot the lattice constant 
a vs. {4>cccc)- Interestingly, in such a "scahng" plot, both classical and quantum results 
are found to collapse nearly on the same straight line, which shows that the lattice constant 
a does not depend on temperature explicitly, but only implicitly, through the temperature 
dependence of the torsional fluctuations. Since the latter dependence is very different in 
classical and in quantum case, the behavior of a is also substantially different. 



In Fig. O, the temperature dependence of the lattice constant h is shown. Similarly to 



the classical case, the error bars for h are larger than those for a. The quantum flattening at 
low temperatures is now less pronounced, and in order to find the limiting zero-temperature 
value of 6, it would be necessary to go to even lower temperatures than 25 K. Here again, 
the quantum curve lies above the classical one at all temperatures. At 25 K, the difference 
between them is about 0.03A, which represents a 0.6 % effect, while at 300 K it decreases 
to about O.O2A. A finite-size effect very similar to the one in the classical case is observed 
here, too. It is visible already at 100 K, where the point for the system with C24 chains 
falls slightly below that for the system with C12 chains, while being strongly pronounced at 
300 K, where the difference is about 0.035A. The same finite-size effect is also visible in the 
"scaling" plot in Fig. n, where moreover at lower temperatures a downward bending of the 
quantum curve can be observed. Comparison to the experimental data Ijl^ in Fig. |13| in 
this case appears to be less good than in the case of a, however, for a detailed comparison 
simulation data for larger system sizes would be required. 

We have also computed the diagonal elastic constants 011,022,033 of the system. Analo- 
gously to the classical case [|ll, we evaluated Cn, C22 from the Parrinello- Rahman fluctuation 
formula ITS, 



Cifc = ^^(eiCfe) A; = l,2,3, (6) 



while for C33 we used the Gusev-Zehnder-Suter formula [|19| in its approximate version suit- 
able for small strain fluctuations 

Cik = -^iPien) {.enCk)'^ , (7) 

n 

where V is the super-cell volume and pi and Ci are the diagonal components of the pressure 
tensor and strain tensor, respectively. This choice of methods was motivated by the finding 
that in the classical case |]l[ significantly smaller statistical errors resulted for C33 from the 
Gusev-Zehnder-Suter formula ([Zp, while for the other elastic constants errors were slightly 
smaller for the Parrinello- Rahman fluctuation formula (|^). Both formulae are classical and 
their use for evaluation of elastic constants also in case of PIMC technique is justified by the 
fact that strain fluctuations are classical objects which have the same values in all Trotter 
slices. The results are shown on Figs, [l^, ^ and |l^. The data for cn exhibit a flattening 
at low temperatures and suggest that the ground state value is reduced with respect to the 
classical one by about 2 GPa. A similar conclusion might be true also for C22, where the 
large statistical error precludes a more detailed comparison. The best results are obtained 
for C33, analogously to the classical case. Here the flattening is clearly seen and the ground 
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state value is reduced due to quantum effects by about 20 GPa. In Fig. [T^, C33 is plotted 
against {(fcccc)- '^^^ quantum points fall close to the line of collapse of the classical points, 
which indicates that a dominant part of the quantum softening of C33 has its origin in the 
finite value of zero-point torsional fluctuations. Not surprisingly, error bars of the quantum 
data are much larger compared to the corresponding classical ones. Our results show that 
for a strongly anisotropic crystal, in might be possible to obtain fairly good results for some 
components of the tensor of elastic constants, while other components might be much more 
difficult to compute. In any case, calculation of elastic constants within the PIMC scheme 
is at present computationally very demanding. 

Finally, we have also studied isotope effects by simulating deuterated PE. In this case, 
we have performed the simulation only at the lowest temperature T = 25 K, with the same 
system size and Trotter number as in case of normal (hydrogenated) polyethylene. The 
results for some quantities are summarized in Table |I[ All three lattice constants are shorter 
in deuterated PE. The largest effect is seen on the lattice constant a, while in case of h its 
relative magnitude is smaller by a factor of 3 and in that of c by a factor of 20. 



IV. CONCLUSIONS 

In this paper, we have demonstrated that for system sizes of several hundred atoms, 
PIMC is a practical method allowing a fully quantum simulation of crystalline systems with 
many different energy scales, like realistic explicit atom models for polymer crystals with 
no constraints on the degrees of freedom. Even in the low-temperature region, where the 
system is close to its ground state, it is possible to calculate lattice constants and internal 
coordinates with a fairly good accuracy. On the other hand, an accurate determination 
of elastic constants is still very difficult. The limitations of the method in a study of the 
thermal expansion of the PE crystal are mainly set by the fact that the finite-size effects in 
the quantum case are more pronounced with respect to the classical one, while at the same 
time it is more difficult to simulate larger systems, because of the extra Trotter dimension. 
It would be very helpful to have for an anisotropic crystal a combined finite-size scaling 
scheme, allowing a simultaneous extrapolation of lattice constants to thermodynamic limit 
and Trotter limit, analogous to the one developed in Ref. [§] for the specific heat of a 
cubic crystal. A prerequisite for such scaling, however, is an availability of high-accuracy 
data. A possible route here might be an improvement on the primitive PIMC algorithm 
by using a better approximation to the density matrix, similar to that used for liquid ^He 
in Ref. 0, allowing a substantial reduction of the Trotter number. In a polymer crystal, 
the stiffest parts of the potential are the bond stretching terms, which formally have a form 
of pair interactions between neighboring atoms. For such pair interactions, it is possible to 
calculate the exact two-body density matrix, either by means of expansion in eigenf unctions, 
or by matrix squaring. This exact form could be tabulated and used in the simulation, while 
the rest of the potential would be treated in a standard way. Such a trick can be expected 
to considerably improve the Trotter convergence, which in turn would enable simulation of 
larger systems and increase the statistical accuracy of the results. 

A detailed comparison to a quasi-harmonic approximation will be done in a forthcoming 
paper pO[. It would also be of interest to perform such an approximation for finite lattices 
in order to clarify the physical origin of the finite-size effects in both classical and quantum 



9 



cases, which still remain to be understood. 

Concerning the physics of the system, we have demonstrated that both in the classical 
and the quantum case, the torsional fluctuations play a central role in the thermal expansion 
of the system. It would be desirable to have an analytical theory of the lateral thermal and 
zero-point expansion, allowing to understand the origin of the anisotropy in both cases. We 
have determined also some local quantities, in particular zero-point fluctuations of internal 
coordinates, like bond lengths and bond and torsional angles, knowledge of which might be 
of relevance for local experimental techniques. 

Finally, we would like to make a few remarks concerning modeling of polymer crystals 
in general, taking into account quantum effects. It seems to us that for this sort of systems, 
force fleld building lags somewhat behind the development of simulation methods, concerning 
in particular the ability to reproduce anharmonic effects correctly. As already pointed out 
in Ref. [|^, in the classical limit, where in principle all phonon modes can contribute to the 
thermal expansion of the system, there is a substantial difference between the properties of 
SLKB force fleld and KDG force fleld as far as the thermal expansion coefficient 



as is concerned. In case of the former one, 0:3 is classically negative at all temperatures, and 
originates dominantly from torsional fluctuations, while for the latter one it vanishes as T — *■ 
0, which points to a large contribution of other modes. Such behavior might perhaps have 
to do with the equilibrium values in the bond stretching and angle bending terms since for 
the KDG force fleld these are substantially different from their average values, which in our 
opinion does not seem to be sufficiently justified. This also demonstrates that the properties 
of a polymer crystal are much more sensitive to the details of the force field than those of a 
liquid. In order to have a systematic control of the important anharmonic properties of a solid 
polymer system, an improvement on the side of the force field building is necessary. Such 
improvement would allow to explore fully the potential of existing classical and quantum 
simulation methods. One first and relatively simple thing to do in order to develop force fields 
suitable for quantum simulations would be to use a quantum quasi-harmonic approximation 
to determine the ground state structure, instead of bare energy minimization, for fitting force 
fields to experimental structures. Obviously, a prerequisite for this is a better experimental 
characterization of the system in the whole range of temperatures, making use of up-to- 
date experimental techniques, like, e.g., x-ray diffraction with synchrotron sources. Such 
techniques should also allow a precise experimental determination of isotope effects on lattice 
constants and thermal expansion, which might in principle help to decide which of the 
abovementioned force flelds provides a better description of the real PE crystal. While for 
the SLKB force fleld the difference between classical and quantum value (the isotope effect 
is also related to this difference) of tends to vanish in the temperature region over 200 
K, where the torsional fluctuations become thermally activated, in case of the KDG force 
fleld a large difference persists up to room temperature f^. Accurate measurements of all 
components of the tensor of elastic constants in a wide region of temperatures would also 
be very interesting and helpful. 
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APPENDIX A 

In this appendix, we present the full form of the force field used in our present quantum 
and previous classical simulation , which is a slightly modified version of the SLKB force 
field [0. While the formal modifications have been described in detail in Ref. here we 
provide the actual numerical values of all parameters. The force field consists of bonded and 
non-bonded interactions. The bonded interactions involve bond stretching, angle bending 
and torsions, as well as off-diagonal, or cross, terms, coupling together the different internal 
coordinates of the chains. The corresponding terms have the form of the following expres- 
sions: 

a) Bond stretching, applying to all C-C and C-H bonds. 

U{r) = h'\r - r,f (8) 

b) Angle bending, applying to all C-C-C, C-C-H and H-C-H angles. 

U{e) = ^k^^^icos 9 - COS ^o)^ (9) 

c) Torsion terms, applying to all C-C-C-C, C-C-C-H and H-C-C-H sequences, 

U{^) = W'''''^{l + cos3^) , (10) 

where ip is the torsional angle {ip = corresponds to cis and (p = n corresponds to trans). 
Throughout the rest of the paper, we use also torsional angle defined as fluctuation of ip 
from the trans value, (p = ip — n. 

d) Bond-angle and bond-bond cross terms, applying to all C-C-C, C-C-H and H-C-H angles, 

U{n, r2, 6) = ki'J{r, - no){cose- cosOo) + K'Jir^ - r2o)(cos0 - cos^o) 

+ ki'J^{n-n^){r2-r2o) , (11) 

where ri and r2 are the bond lengths of the I- J and J-K bonds adjacent to an UK bond 
angle 6. 

e) One center angle-angle cross terms, applying to all pairs of bond angles about a tetrahedral 
carbon atom sharing a common bond, 

U{ei,e2) = G^-^'-^^ {cos 01 - cos^t)(cos^2 - cosOt) , (12) 

where 61,62 are the JIK and JIL bond angles, respectively, and I is the tetrahedral carbon. 
Here, 6t is the tetrahedral angle. 

f) Two-center angle-angle cross terms, applying to all C-C-C-C, C-C-C-H and H-C-C-H 
sequences, 

U{ip, 61, 62) = F^'-^^'-^ cosip{cos6i - cos6t)icos62 - cos6t) , (13) 

where (p is the torsional angle corresponding to the sequence IJKL, 61,62 are the UK and 
JKL bond angles and 6t is again the tetrahedral angle. Numerical values of all the parameters 
in the above expressions are contained in the table Tab.||. 
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The non-bonded interaction has the form U{r) = Ae~^^ — Cr~^ between atoms on dif- 
ferent chains and atoms on the same chain separated by more than two bonds (1-2 and 
1-3 interactions are excluded). Numerical values of all the parameters for all three pairs 
of atoms (C-C, H-H and C-H) are contained in the table Tab.|T|. Details concerning cutoff 
and long-range corrections can be found in Ref. ||1|. 

APPENDIX B 

In this appendix we discuss in detail the origin of the quantitative discrepancy between 
the slope of the c vs. |(((/>o ~ line as obtained from formula and from our quantum 
simulation data (Fig. p^ ). To this end, we have to look more closely to the way the formula 
has been derived in Ref. |0| . The model used assumes a single chain consisting of rigid C-C 
bonds with length a and rigid C-C-C angles Occc = n — a, subject to torsional deformations 
only. The expression (^ then evaluates c simply as the average end-to-end distance between 
carbon atoms separated by four bonds. Let us denote by r, r the instantaneous positions 
of the two atoms, and by tq, their equilibrium positions around which they oscillate with 
instantaneous fluctuations A, A'. The squared end-to-end distance is then given by 

(f - = (fo - ro + A - A'Y = (fo - To)' + (A - A? , (W) 

where the cross term has vanished, because the fluctuations A, A' are orthogonal to tq " 
Tq (since the bond lengths and angles are assumed to be rigid, the only possible small 
displacements of carbon atoms are those perpendicular to the plane of the unperturbed 
all-trans chain). This expression contains apart from the squared distance between the 
equilibrium positions (ro — Tq)^, which is directly related to the lattice constant c, also a 
fluctuation term. Neglecting this term for a relatively short segment of chain, as done in 
Ref. 1^, results in underestimating of the thermal contraction. In order to improve on 
this, one has to consider a longer segment of the chain, which would, however, require 
a knowledge of correlation functions between torsions displaced by more than one bond. 
This is clearly impossible within the simple single-chain model assumed in Ref. ||^, where 
the crystal environment of the chain is neglected entirely. This together with a separable 
torsional potential classically leads to vanishing of all correlation functions between displaced 
torsions, and in turn to coiling of longer segments of a chain. In a crystal it is just the external 
non-bonded field of all other chains acting directly on the absolute coordinates of the atoms 
of a given chain (rather than on the torsional angles, which by their very nature have a 
character of relative coordinates), which induces non-zero correlations between the torsions 
in the chain. These correlations reflect the existence of the underlying 3D crystal structure 
with its translational long-range order and result in an overall coherent shortening of the 
chain, instead of its coiling. The quantitative agreement of the formula with experiment, 
as stated in Ref. P], thus appears to be accidental and arises due to compensation of two 
effects: neglecting the crystal field contribution lowers the effective torsional constant by 
a factor of two which in turn increases the torsional fluctuations and compensates for the 
lower value of the proportionality constant in the expression (|^). 

Now we derive a few formulas similar to (^, taking into account progressively longer 
segments of a chain. In general, we are interested in calculating the end-to-end distance of a 
segment of a chain consisting of = 2m bonds. The corresponding vector can be expressed 
as 
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Rn = J2I[T, (a, 0,0)^ 



k=l i=l 

where the 3x3 matrices T,- are defined as follows 



(15) 



cos a sin a 

sin a cos (pi — cos a cos 0j sin (p. 
sin a sin — cos a sin — cos ( 



(16) 



angles (pi being the torsional angles. Evaluating {\Rn\) and keeping just terms up to second 
order in the angles (pi we find the following approximations c'-™'-' = ^(|-R2m|) to the contracted 
lattice constant c (the actual calculation has been performed by Mathematica) 



.(2) 



.(3) 



.(4) 



.(5) 



.(6) 



2(l^4|) 



^(l^sl) 

(0005)) 

^(l^iol) 



Co 
Co 
Co 

= Co 



isin2|((0o^)-(0o0i)) 
i sin^ |(4(0o') - HMi) + 2(0002) - (0003)) 



a 



16 



sm 



(10(00^) - 14(0o0l) + 8(0002) - 5(0003) + 2(0004) 



1 - ^sin^ |(2O(0o') - 3O(0o0i) + 2O(0o02) - 14(0o03) + 8(0o04) 



5(0005) +2(0006) - (0007)) 



6(l^i2l) 



Co 



1 a 

1 - — sin^ -(35(0o') - 55(0o0i) + 4O(0o02) - 3O(0o03) + 2O(0o04) 
36 / 



- 14(0005) + 8(0006) - 5(0007) + 2(0008) " (0009)) 



(17) 



where Cq = 2a cos ^ is the lattice constant of the unperturbed chain. 

We introduce now off-plane displacements Zi of carbon atoms along local ^-axes defined 
for each C atom by the unit vector ki = {Si x aj+i)/ sin a, where Si is the vector connecting 
carbon atoms i — 1 and i. Vectors ki are perpendicular to the plane of the all-trans chain and 
their directions alternate below and above the plane. It can then be easily shown |^ that the 
torsional angles (pi defined via cos0j = — (aj_i x Oj) ■ (cij x aj_|_i)/a^ sin^ a, or, alternatively, 
sin0j = aj_i • {di X aj_|_i)/a^ sin^ a, can be in first order in displacements Zi expressed as 
(pi = (— ^i_2 — + Zi + Zi+i)/asma. Upon substitution of this latter expression to (|17D, 
we find 



^(m) 



1 {izo-Z2mf) 

2rn? Co 



where 



Co 1 



((^0-^2)^)' 



2cl 



(18) 



(19) 



is the limiting m oo value of the contracted lattice constant c. The last result has 
been derived also in Ref. by means of a simple geometrical argument, expressing c as a 
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projection of the segment consisting of two C-C bonds, whose end-to-end length is cq, on 
the plane of the unperturbed chain. 

Applying the formula for c^^-* corresponding to a segment consisting of six C-C bonds to 
our simulation data, we find a slope which is still larger by about 50 % than the theoretical 
one. A characteristic property of all the expressions (p!7| ) is that the coefficients of torsion 
correlation functions increase in magnitude with the length of the segment. The correct lim- 
iting result for thermal contraction thus emerges progressively as a consequence of a delicate 
cancellation between the terms, which means that a knowledge of correlation functions at 
many different displacements with very good accuracy would be required. Such behavior 
reflects the lack of existence of a preferred crystal direction for a chain in the formulation 
employing exclusively torsional angles. In contrast to this, upon substitution of the true off- 
plane displacements of the atoms into these expressions, the limiting exact result emerges in 
a very simple form, which has a clear geometrical interpretation and contains just a single 
correlation function between second neighbors. This demonstrates that while the relation 
between the torsional fluctuations and the lattice constant c provides a very useful insight 
for qualitative understanding of the thermal contraction of c, it would be at the same time 
hard to push these considerations to a truly quantitative level. 
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FIGURES 



FIG. 1. Temperature dependence of the average fluctuation ^ {{5rcc)'^) of the C-C bond length, 
in classical and quantum case, shown for different system sizes. In this and most of the following 
figures, statistical error bars are shown, lines are for visual help only. In all figures, the same symbol 
is used for quantum results corresponding to different values of the Trotter number P. When at a 
given temperature the results for different values of P are indistinguishable within the statistical 
error, like in this figure, the Trotter numbers are not indicated explicitly. 

FIG. 2. Temperature dependence of the average fluctuation {{SrcnY) of the C-H bond 
length, in classical and quantum case, shown for different system sizes. When the quantum results 
at a given temperature exhibit a pronounced dependence on the Trotter number P, like in this 
figure, the corresponding values of P are indicated next to the symbols. 



FIG. 3. Temperature dependence of the average fluctuation \/ {{^OcccY) of the C-C-C bond 
angle, in classical and quantum case, shown for different system sizes. 



FIG. 4. Temperature dependence of the average fluctuation a/ {{^OhchY) of the H-C-H bond 
angle, in classical and quantum case, shown for different system sizes. 



FIG. 5. Temperature dependence of the average torsional angle fluctuation y {4>'cccc) ^om 
the trans minimum, in classical and quantum case, shown for different system sizes. 

FIG. 6. Temperature dependence of the internal energy per unit cell of quantum system with 
Ci2 chains, shown for different values of the Trotter number P, together with an extrapolation to 
P ^ oo. Note the strong dependence on the Trotter number P. 

FIG. 7. Temperature dependence of the lattice parameter c, in classical and quantum case. 



shown for different system sizes, together with the experimental data [17|. 



FIG. 8. Lattice parameter c vs. {(pcccc)^ classical and quantum case, shown for different 
system sizes. Note that also in the quantum case, the dependence is almost linear over the whole 
temperature range. The values of temperature, which is a parameter of the plot, are indicated 
next to the symbols. 

FIG. 9. Temperature dependence of the normalized correlation function of fluctuations 

of two neighboring torsional angles from the trans minima, in the classical and quantum case, for 
the system with C12 chains. Note the pronounced dependence on the Trotter number P. 
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FIG. 10. Lattice parameter c vs. ^{{(/)o — (f>i)'^) in the quantum case for system with C12 chains. 
The values of temperature, which is a parameter of the plot, are indicated next to the symbols. 
The points can be fitted by a line c = 2.5327 - 1.275 x 10-'^i(((/)o - cl)i)^). 

FIG. 11. Temperature dependence of the lattice parameter a, in classical and quantum case, 
shown for different system sizes, together with the experimental data [17|. 

FIG. 12. Lattice parameter a vs. {(pcccc)^ classical and quantum case, shown for different 
system sizes. The values of temperature, which is a parameter of the plot, are indicated next to 
the symbols. Note the collapse of both classical and quantum results on almost the same straight 
line over the whole temperature range. 

FIG. 13. Temperature dependence of the lattice parameter b, in classical and quantum case. 



shown for different system sizes, together with the experimental data [17|. Note the strong finite-size 
effect at higher temperatures for both classical and quantum results. 

FIG. 14. Lattice parameter b vs. {(fcccc)^ ™ classical and quantum case, shown for different 
system sizes. The values of temperature, which is a parameter of the plot, are indicated next to 
the symbols. Note the strong finite-size effect at higher temperatures for both sets of data as well 
as the downward bending of the quantum curve at low temperatures. 

FIG. 15. Elastic constant cn as a function of temperature, in both classical and quantum case, 
shown for different system sizes. 

FIG. 16. Elastic constant C22 as a function of temperature, in both classical and quantum case, 
shown for different system sizes. 

FIG. 17. Elastic constant C33 as a function of temperature, in both classical and quantum case, 
shown for different system sizes. 

FIG. 18. Elastic constant C33 vs. {4>cccc)^ ™ both classical and quantum case, for system with 
C12 chains. The values of temperature, which is a parameter of the plot, are indicated next to the 
symbols. Note that the quantum results fall close to the line passing through the classical results. 
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TABLES 





normal PE 


deuterated PE 


difference 


a 


s 

7.201 ± 0.003^ 


5: — 

7.170 ±0.0031 


-0.43 % 


b 


4.928 ± 0.0021 


4.921 ± 0.0021 


-0.14 % 


c 


2.5297 ± O.OOOlA 


2.5292 ± O.OOOll 


-0.02 % 


V {4>cccc) 


5.59 ±0.01° 


5.26 ± 0.015° 


-5.9 % 


VmHCHY) 


8.28 ±0.001° 


7.17 ±0.003° 


-13.4 % 




0.074 ± O.OOOOll 


0.065 ± O.OOOOll 


-12.1 % 



TABLE 1. Isotope effect: deuterated PE compared to normal (hydrogenated) PE. All values 
correspond to the system with C12 chains and P = 144 at T = 25 K. No extrapolation to P ^ 00 
has been performed here. 





1.53 


fc^^^ [kcal/(moll)] 





k^^ [kcal/(moll2)] 


617.058 


k^g^ [kcal/(moll2)] 







1.09 


f^rhe^ [kcal/(moll)] 


-34.0818 


k^^ [kcal/(moll2)] 


654.455 


kg^g^ [kcal/(moll)] 





^HCH [kcal/mol] 


76.952 


fc.^^f [kcal/(mol 12)] 







107.899 


[kcal/(moll)] 


-54.494 


f^cCH [kcal/mol] 


85.726 


fc^^^' [kcal/(moll2)] 


25.9337 




109.469 


qCC:CH [kcal/mol] 


-6.0034 


fc^^^ [kcal/mol] 


107.446 


qCC-.hh [kcal/mol] 


-3.6409 




110.999 


qCH-.cc [kcal/mol] 


2.9127 


yHCCH [kcal/mol] 


0.2776 


qCH-.ch [kcal/mol] 





yCCCH [kcal/mol] 


0.2776 


pH:CC:H [kcal/mol] 


-16.0398 


yCCCC [kcal/mol] 


0.2776 


pC:CC:H [kcal/mol] 


-15.5343 






pC:CC:C [kcal/mol] 


-33.3664 



TABLE IL Parameters of the bonded interaction of the force field. The left column contains 
the coefficients of terms diagonal in the internal coordinates of the chains, while the right one 
contains the coefficients of the off-diagonal terms. 



atoms 


A [kcal/mol] 


B [1-1] 


C [A^' kcal/mol] 


C-C 


14889.0 


3.089 


639.58 


H-H 


2640.2 


3.739 


27.39 


C-H 


4300.9 


3.416 


137.44 



TABLE in. Parameters of the non-bonded interaction of the force field. 
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